Excitation spectrum and high energy plasmons in single- and multi-layer graphene 
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In this paper we study the excitation spectrum of single- and multi-layer graphene beyond the 
Dirac cone approximation. The dynamical polarizability of graphene is computed using a full n- 
band tight-binding model, considering the possibility of inter-layer hopping in the calculation. The 
effect of electron-electron interaction is considered within the random phase approximation. We 
further discuss the effect of disorder in the spectrum, which leads to a smearing of the absorption 
peaks. Our results show a redshift of the 7r-plasmon dispersion of single-layer graphene with respect 
to graphite, in agreement with experimental results. The inclusion of inter-layer hopping in the 
kinetic Hamiltonian of multi-layer graphene is found to be very important to properly capture the 
low energy region of the excitation spectrum. 
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I. INTRODUCTION 

One of the main issues in the understanding of the 
physics of graphene is the role played by electron-electron 
interaction^ Several collective modes as low and high en- 
ergy plasmons, as well as plasmarons, are a consequence 
of electronic correlations and have been measured in this 
material. The high energy 7r-plasmons have been ob- 
served in electron energy-loss spectroscopy (EELS)jSr— 
inelastic X-ray scattering (IXS)£ or optical conductivity^ 
Recently, a plasmaron mode (which is a result of coupling 
between electrons and plasmons) has been measured in 
angle- resoved photoemission spectroscopy (ARPES).— 

At low energies, long range Coulomb interaction leads, 
in doped graphene, to a gapless plasmon mode which 
disperses as uj p i ~ ^fq^ and which can be described 
theoretically within the random phase approximation 
(RPA))2~— The low energy linear dispersion relation of 
graphene is at the origin of a new series of collective 
modes predicted for this material and which do not 
exist for other two-dimensional electron gases (2DEG), 
as inter- valley plasmons^ or linear magneto-plasmonsr^ 
which can be described within the RPA as well. For un- 
doped graphene, the inclusion of ladder diagrams in the 
polarization can lead to a new class of collective modest 
as well as to an excitonic instability^^— 

However, much less is known about the high energy 
7r-plasmon which, in the long wavelength limit, has an 
energy of the order of 5-6 eV, and which is due to the 
presence of Van Hove singularities in the band disper- 
sion. For single- layer graphene (SLG), this mode has 
been studied by Stauber et &\M and by Hill et ali24 in 
the RPA. Yang et al. have included excitonic effects and 
found a redshift of the absorption peak^ leading to a 
better agreement with the experimental results^ Here 
we extend those previous works and study the excitation 
spectrum of SLG and multi-layer graphene (MLG) from 
a tight-binding model on a honeycomb lattice. By means 
of the Kubo formula, the non-interacting polarization 
function n(q, uj) is obtained from the numerical solution 
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FIG. 1. Atomic structure of ABA- and ABC-stacked multi- 
layer graphene. The intra-layer (t) and inter-layer (71 and 73) 
hopping amplitudes are considered, as explained in the text. 



of the time-dependent Schrddinger equation. Coulomb 
interactions are considered in the RPA, the validity of 
which is discussed. We also consider the effect of disor- 
der in the system, which lead to a considerable smearing 
of the Van Hove singularities in the spectrum. Our re- 
sults show a redshift of the 7r-plasmon mode in graphene 
with respect to graphite, as it has been observed in the 
experiments!^ Furthermore, the inclusion of inter-layer 
hopping is found to be very important to capture the low 
energy region of the spectrum in MLG. 



The paper is organized as follows. In Sec. |H] we de- 
scribe the method that we use to compute the dynamical 
polarization function of SLG and MLG. In Sec. Mil we 
give results for the excitation spectrum of SLG, consid- 
ering the effect of disorder and electron-electron interac- 
tion. The spectrum of MLG is described in Sec. IIVI In 
Sec. |V] we compare our results to recent experimental 
data. Our main conclusions arc summarized in Sec. I VII 
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II. DESCRIPTION OF THE METHOD 



The tight-binding Hamiltonian of a MLG is given by 

(1) 



JVlayer N laycr -1 

H = E H t+ E K 

1=1 i=i 



where Hi is the Hamiltonian of the Vth layer of graphene, 

H i = - E (*Mi a /,i 6 w + h - c ) + E w m4; c m> ( 2 ) 

where a| 4 creates (annihilates) an electron on sub- 
lattice A (B) of the i'th layer, and ti ij is the nearest 
neighbor hopping parameter. The second term of Hi ac- 
counts for the effect of an on-site potential vi.i, where 
nij = CtjCii is the occupation number operator. In the 
second term of the Hamiltonian Eq. (JTJ) , H[ describes the 
hopping of electrons between layers I and I + 1. In nature 
there are two known forms of stable stacking sequence in 
bulk graphite, namely ABA (Bernal) and ABC (rhom- 
bohedral) stacking, and they are schematically shown in 
Fig. [TJ For a MLG with an ABA stacking, H[ is given 
by 

H[ = - 7l [<4,j b i+i,j + hx - E [ 6 L°*+u' + h - c - 

3 3,3' 

(3) 

where the inter-layer hopping terms 71 and 73 are shown 
in Fig. [1] Thus, all the even layers (7 + 1) are rotated with 
respect to the odd layers (I) by +120°. The difference be- 
tween ABA and ABC stacking is that, the third layer(s) 
is rotated with respect to the second layer by —120° 
(then it will be exactly under the first layer) in ABA 
stacking, but by +120° in ABC stacking^ In this paper 
we use the hopping amplitudes t = 3 eV, 71 = 0.4 eV 
and 73 = 0.3 eV^i The spin degree of freedom con- 
tributes only through a degeneracy factor and is omitted 
for simplicity in Eq. ([I]). In our numerical calculations, 
we use periodic boundary conditions in the plane (XY) 
of graphene layers, and open boundary conditions in the 
stacking direction (Z). 

The dynamical polarization can be obtained from the 
Kubo formula^ as 

U(d,uj) = v J o dTe WT ([p(q,r),p(-q,0)]) ) (4) 

where V denotes the volume (or area in 2D) of the unit 
cell, p (q) is the density operator given by 

N layer 

p(i)= E E c m c m ex p (*q • r u) > ( 5 ) 

1=1 i 

and the average is taken over the canonical ensemble. For 
the case of the single-particle Hamiltonian, Eq. (0} can 
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FIG. 2. 2D Brillouin zone of SLG. For undoped graphene, the 
valence and conduction bands touch each other at the vertices 
of the hexagon, the so called Dirac points (K and K'). The 
Van Hove singularity lies at the M point, and we have defined 
9 as the angle between the wave-vector q and the fc x -axis. 



be written as^ 

2 f°° 

H(q, W ) = --y dre^ 

xlm (<p\ n F (H) e lHT p (q) e"^ [1 - n F (H)\ p (-q) \<p) , 

(6) 

where np (H) = e/HH is the Fermi-Dirac distribu- 
tion operator, /? = 1/ksT where T is the temperature 
and ks is the Boltzmann constant, and p, is the chemi- 
cal potential. In the numerical simulations, we use units 
such that h = 1, and the average in Eq. ^ is performed 
over a random phase superposition of all the basis states 
in the real space, i.e. ) 29 i 30 



|v)=E a M c Mi°) 



(7) 



where an are random complex numbers normalized as 

i \ a i,i\ 2 = 1- By introducing the time evolution of two 
wave functions 

|^ (q,r)) = e~ iH T [1 - n F (H)] p (-q) \<p) , (8) 
\<P2 (r)> = e~ iHT n F (H) \<p) , (9) 

we get the real and imaginary part of the dynamical po- 
larization as 

2 f°° 

Ren(q,w) = -— J dTco8(uT)Jm{(p2{r)\p(q)\(pi(T)), 



Imn (q, uj) = -— J dr sin(wr) Im (ip 2 (r) \p (q)| ipi (r)) , 

(10) 

The Fermi-Dirac distribution operator np (H) and the 
time evolution operator e~ iHr can be obtained by the 
standard Chcbyshcv polynomial decomposition^ 
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For the case of SLG, we will further compare the polar- 
ization function obtained from the Kubo formula Eq. Q , 
to the one obtained from the usual Lindhard function^ 
Notice that this method can be used to calculate the 
optical conductivity of graphene beyond the Dirac cone 
approximation i 29 i 32 For pristine graphene, the dynamical 
polarization obtained from the Lindhard function using 
the full 7r-band tight-binding model is£ ' 23 ' 24 



n(q,w) 



d 2 k 



E /-'(k,q) 



(2tt) z Jbz 

np [E s (k)] — np E s ' (k + q) 



E s (k) - E s ' (k + q) + u + iS ' 



(11) 

where the integral is over the Brillouin zone, g s = 2 is 
the spin degeneracy, E^ (k) = ±t \(j>k\ — \i is the energy 
dispersion with respect to the chemical potential, where 



1 + 2e i3fcW2 cos f^ fe 



(12) 



a = 1.42A being the in-plane carbon-carbon distance, 
and the overlap between the wave-functions of the elec- 
tron and the hole is given by 



/±(k,q) 



l±Rc 



I0k+q| 



(13) 



In the RPA, the response function of SLG due to 
electron-electron interactions can be calculated as 



x(q,w) = 



n(q,w) 



V(</)n(q,w)' 



where V (q) = is the Fourier component of the 

Coulomb interaction in two dimensions, in terms of the 
background dielectric constant k, and 



e(q,w) = l-V(q)U{q,u) 



(15) 



is the dielectric function of the system. We will be in- 
terested on the collective modes of the system, which 
are defined from the zeroes of the dielectric function 
[e(q, uj) = 0]. The dispersion relation of the collective 
modes is defined from 

Re e(q, oj p i) = 1 - V(g)II(q, u: pl ) = 0, (16) 

which leads to poles in the response function ([14")) . The 
damping 7 of the mode is proportional to Im n(q, u> p i), 
and it is given by 

Imn(q,^) 
7 fRcn(q,w)| ' V 

For MLG, the response function is calculated as (we 
use q z 



X3D (q,w) 



IL 3D (q, w) 



V(q)F(q)Il 3D (q,w)d' 



(18) 



where d = 3.35A is the inter-layer separation. Because 
we use open boundary conditions in the stacking direc- 
tion, we define the form factor F (q) as 



F(q) = 



1 



layer 



E * 

l,i' =1 



q\l-V\d 



(19) 



The expression (fT5]) assumes that the polarization of each 
layer is the same, and it is exact in two different limits: 
bilayer graphene and graphite. Notice that a similar ef- 
fective form factor has been used to study the loss func- 
tion of multiwall carbon nanotubes^ 3 - Eq. (fH)]) coincides 
with the commonly used form factor for a multi-layer 
system with an infinite number of layers^ 



F(q) 



N,, 



E' 



-q\l-l'\d 



(20) 



where in this last case the periodicity ensures that F (q) 
is independent of layer index I, with the asymptotic be- 
havior F(q) = sinh(grf)/[cosh(grf) — 

A crucial issue is the value of the dielectric constant n 
for each of the cases considered, because it encodes the 
screening due to high energy (a) bands which are not ex- 
plicitly considered in our calculation. A good estimation 
for it can be obtained from the expression 3 ^ 



K(q) 



«i 



+ l-(«i- l)e-« L 



Ki + 1 + (ki - l)e-i L 



Ki, 



(21) 



where k\ w 2.4 is the dielectric constant of graphite, L = 
d m + (Ni a y er — 1) d is the total height of the multi- layer 

(14) 

system in terms of the number of layers Ni ayer and the 
height of a monolayer graphene d m 2.8 A. As expected, 
Eq. ([2~1"T) gives k = 1 for SLG at q — >• and k = «i for 
graphite. 

We notice that the accuracy of the numerical results 
for the polarization function Eq. (|10|) is mainly deter- 
mined by three factors: the time interval of the propa- 
gation, the total number of time steps, and the size of 
the sample. The maximum time interval of the propaga- 
tion in the time evolution operator is determined by the 
Nyquist sampling theorem. This implies that employing 
a sampling interval At = 7r/maxi \Ei\, where are the 
cigenencrgies, is sufficient to cover the full range of en- 
ergy eigenvalues. On the other hand, the accuracy of the 
energy eigenvalues is determined by the total number of 
the propagation time steps (N T ) that is the number of 
the data items used in the fast Fourier transform (FFT). 
Eigenvalues that differ less than AE = it/N t At cannot 
be identified properly. However, since AE is proportional 
to N^ 1 we only have to double the length of the calcu- 
lation to increase the accuracy by the same factor. The 
statistic error of our numerical method is inversely pro- 
portional to the dimension of the Hilbert spaced and 
in our case (the single particle representation), it is the 
number of sites in the sample. A sample with more sites 
in the real space will have more random coefficients (a/^) 
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FIG. 3. (Color online) — Im Il(q, u) for SLG in the clean 
limit for different values of wave- vector q. Plots (a), (c) and 
(e) correspond to a wave- vector q along F-K, whereas (b), 
(d) and (f) are of a q parallel to the T-M direction. The 
angle 9 is defined in Fig. [2] In the numerical integration of 
Lindhard function in Eq. pip , we use 2 x 10 8 Monte Carlo 
points (k) in the first Brilion. The sample size of SLG used 
in the numerical calculation of Kubo formula in Eq. ([6]) is 
4096 x 4096. 



in the initial state | ip) , providing a better statistical repre- 
sentation of the superposition of all energy eigenstates<^ 
Similar algorithm has been successfully used in the 
numerical calculation of the electronic structure and 
transport properties of single- and multi-layer graphene, 
such as the density of states (DOS), or dc and ac 
conductivities! 29 ' 36 ' 37 The main advantage of our algo- 
rithm is that different kinds of disorders and boundary 
conditions can be easily introduced in the Hamiltonian, 
and the computer memory and CPU time is linearly pro- 
portional to the size of the sample, which allows us to do 
the calculations on a sample containing tens of million 
sites. 



III. EXCITATION SPECTRUM OF 
SINGLE-LAYER GRAPHENE 

The particle-hole excitation spectrum is the region 
of the energy-momentum space which is available for 
particle-hole excitations. For non-interacting electrons, 
it is defined as the region where Im II(q, cj), as given 
by Eq. (0| or (fTTj) . is non-zero<2I The linear low energy 




FIG. 4. (Color online) Density of states for SLG considering 
different kind of disorder. The left inset shows a zoom of 
the DOS near the Dirac point (E = 0), whereas the right 
hand side inset shows the disorder broadening of the Van 
Hove singularity at E = t. The numerical method used in the 
calculation of DOS is discribed in Ref. |2SL and the sample 
size of SLG is 4096 x 4096. 



dispersion relation of graphene as well as the possibility 
for inter-band transitions lead to a rather peculiar ex- 
citation spectrum for SLG as compared to the one of a 
two-dimensional electron gas (2DEG) with a parabolic 
band dispersion^ Here we focus on undoped graphene 
(fx = 0), for which only inter-band transitions are al- 
lowed. In Fig. [3] we plot Im II(q, to) for different wave- 
vectors at T = 300/^ (which is the temperature that we 
will use from here on in our results). The first thing 
one observes is the good agreement between the results 
obtained from the Kubo formula Eq. ((4]), as compared 
to Lindhard function Eq. what proofs the valid- 

ity of our numerical method. Furthermore, for the small 
wave- vector used in Fig. G3a)-(b), the results are well de- 
scribed by the Dirac cone approximation , 10 ' 1 1 but only 
at low energies, around uj ~ vpq, where v-p = Sat/2 is 
the Fermi velocity near the Dirac points. In particular, 
the continuum approximation cannot capture the peaks 
of ImII(q, u>) around tu m 2t. These peaks are related 
to particle-hole excitations between states of the valence 
band with energy E w —t and states of the conduction 
band with energy E sa t, which contribute to the polar- 
ization with a strong spectral weight due to the enhanced 
density of states at the Van Hove singularities of the ir- 
bands (see Fig. 

Second, for larger wave-vectors [Figs. EKc)-(f)] one 
observes strong differences in the spectrum depending 
on the orientation of q, effect which has been discussed 
previousl y 23 ' 39 If q is along the T-K direction, there is 
a splitting of the peak associated to the Van Hove sin- 
gularity at lo ~ 2t. At low energies, we also observe a 
finite contribution to the spectral weight to the left of the 
lu vpq peak for momenta along the T-M direction [plots 
Figs. EJd) and (f)]. Finally, trigonal warping effects are 
important as we increase the magnitude of |q|, due to 
the deviation of the band dispersion with respect to the 
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FIG. 5. (Color online) — Im n(q,aj) for different kinds of disorder and for different values and orientation of the wave-vector q. 
In all the plots, the results using the Kubo formula Eq. Q are compared to the Dirac cone approximation. The sample size 
of SLG is 4096 x 4096. 



linear cone approximation. As a consequence, the con- 
stant energy contours are not any more circles around the 
Dirac points, but present a triangular shape. The con- 
sideration of this effect leads to a redshift of the w~«f? 
peak with respect to the Dirac cone approximation, as 
seen clearly in Fig. (He). 

Once we have discussed the clean case, we consider the 
effect of disorder on the excitation spectrum as explained 
in Sec. HU Two different kinds of disorder are consid- 
ered: random local change of on-site potentials and ran- 
dom renormalization of the hopping, which correspond 
to the diagonal and off-digonal disorders in the single- 
layer Hamiltonian Eq. ([2]), respectively. The former acts 
as a chemical potential shift for the Dirac fermions, i.e., 
shifts locally the Dirac point, and the later rises from the 
changes of distance or angles between the p z orbitals. In 
Fig. @]wc show the DOS of SLG for different kinds and 
magnitudes of disorder. The DOS for clean graphene has 
been plotted by using the analytical expression given in 
Ref. |4fJ. The DOS of the disordered systems are calcu- 
lated by Fourier transform of the time-dependent corre- 
lation functions^ 



P(e) 



1 

2^ 



(<4 



iHI 



\f)dt, 



(22) 



with the same initial state \ip) defined in Eq. 0. As 



shown in Ref. [29j, the result calculated from a SLG with 
4096 x 4096 lattice sites matches very well with the ana- 
lytical expression, and here we use the same sample size 
in the disordered systems. We consider that the on-site 
potential Vi is random and uniformly distributed (inde- 
pendently on each site i ) between —v r and +v r . Sim- 
ilarly, the in-plane nearest-neighbor hopping tij is ran- 
dom and uniformly distributed (independently on sites 
i,j) between t — t r and t + t r . The main effect is a smear- 
ing of the Van Hove singularity at E — t, as observed in 
the right hand side inset of Fig. |U 

The effect of disorder is also appreciable in the non- 
interacting excitation spectrum of the system, as shown 
by Fig. [5J A broadening of the w « upg and uj w 2t 
peaks is observed in all the cases. Furthermore, disorder 
leads to a slight but appreciable redshift of the peaks with 
respect to the clean limit. This effect is more important 
for higher wave- vectors, as it can be seen in Fig. E][c)-(d). 
Finally, the disorder broadening of the peaks leads in all 
the cases to a transfer of spectral weight to low energies 
(below lo = v-pq), as it is appreciable in Fig. E^a)-(d). 

The next step is to consider both, disorder and 
electron-electron interaction in the system. In the RPA, 
the response function is calculated as in Eq. (|14p . The 
results are shown in Fig. [6j where — Im x(q, uj) is plot- 
ted for the same wave- vectors and disorder used in Fig. 
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FIG. 6. (Color online) — Im x(*Ji CJ ) f° r the same values of q and disorder as in Fig. 




co(eV) 
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FIG. 7. (Color online) Re e(q, u) for the same values of q and 
disorder as in Fig. [5] and [5] 



We observe that the Dirac cone approximation (black 
line) captures well the low energy region of the spectrum. 
However, the large peak at oj ~ 2t cannot be captured by 
the continuum approximation. They are due to a plas- 
mon mode associated to transitions between electrons in 
the saddle points of the 7r-bands. Strictly speaking, those 
modes cannot be considered as fully coherent collective 



modes, as for example, the low energy y^-plasmon which 
is present in doped graphene^ For doped graphene, the 
acoustic y^-plasmon is undamped above the threshold 
u = v-pq until it enters the inter-band particle-hole con- 
tinuum, when it starts to be damped and decays into 
electron-hole pairs. However, the 7r-plasmon, although it 
corresponds to a zero of the dielectric function as it can 
be seen in Fig. [3 it is a mode which lies inside the con- 
tinuum of particle- hole excitations: — Im n(q, uj p {) > at 
the 7r-plasmon energy ui p i , and the mode will be damped 
even at q — > 0. In any case, it is a well defined mode 
which has been measured experimentally for SLG and 
MLGi£~— Coming back to our results, notice that the 
height of the peaks is reduced when the effect of disorder 
is considered, although the position is unaffected by it. 
For small wave-vectors, this mode is highly damped due 
to the strong spectral weight of the particle-hole excita- 
tion spectrum at this energy, as seen e.g. by the peak 
of -Im II(q,w) at uj = 2t in Fig. Eta)-(b). The posi- 
tion of the collective modes can be alternatively seen by 
the zeroes of the dielectric function Eq. (fll)]) . which is 
shown in Fig. [7] Notice that the Dirac cone approxima- 
tion (solid black lines in Fig. [7]) is completely insufficient 
to capture this high energy 7r-plasmon, which predicts 
always a finite Re e(q, ui). As for the polarization, we 
see that disorder lead to an important smearing of the 
singularities of the dielectric function, as seen in Fig. [7] 
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FIG. 8. (Color online) (a,b) Density of states of ABA- (left 
panels) and ABC-stacked (right panels) multilayer graphene. 
A zoom of the DOS around the Dirac point (E = 0) is shown 
in (c,d), and around the Van Hove singularity (E = t) is 
shown in (e,f). The sample sizes of each layer in MLG are: 
4096 x 4096 atoms in bilayer; 3200 x 3200 in trilayer; 2048 x 
2048 in 5 layers; 1600 x 1600 in 10 layers and 800 x 800 in 50 
layers. 



Finally, we mention that the application of our method to 
even higher wave- vectors and energies as the ones consid- 
ered in the present work, should be accompanied by the 
inclusion of local field effects (LFE) in the dielectric func- 
tion, which are related to the periodicity of the crystalline 
latticed In fact, for SLG and for wave- vectors along the 
zone boundary between the M and the K points of the 
Brillouin zone (see Fig. [2]), the inclusion of LFEs leads to 
a new optical plasmon mode at an energy of 20-25 eV42 



IV. EXCITATION SPECTRUM FOR 
MULTI-LAYER GRAPHENE 

In the following, we study the excitation spectrum and 
collective modes of MLG. For this, we consider not only 
the Coulomb interaction between electrons on different 
layers, but also the possibility for the carriers to tun- 
nel between neighboring layers, as described in Sec. [TXJ 
The importance of considering inter-layer hopping has 
been already shown in the study of screening properties 
of MLG42 First, we see that the results are sensitive to 
the relative orientation between layers. In Fig. |S]we show 
the density of states for ABA- and ABC- stacked MLG 
(see Fig. [1] for details on the difference between those 
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FIG. 9. (Color online) Band structure of ABA- and ABC- 
stacked trilayer graphene. Left panel: the red dashed lines 
indicate the position of the jump in DOS of ABA-stacked 
trilayer graphene at \E\ « 0.19t. Right panel: the red dashed 
lines indicate position of the peaks in DOS of ABC-stacked 
trilayer graphene at E = and \E\ ~ 0.12f. 



two orientations). As seen in Fig. rH^c)-(d), all the MLGs 
present a finite DOS at E = 0, contrary to SLG which 
has a vanishing DOS at the Dirac point. The main differ- 
ence between the two kinds of stacking is that for ABC 
there is a central peak together with a series of satellite 
peaks around E = [Fig. Efd)], whereas for ABA the 
DOS follows closer the behavior of the SLG [FigEJc)]. 
The different structure in the DOS can be understood by 
looking at Fig. |H1 where we show the low energy band 
structure of a trilayer graphene with ABA [Fig. Ela)] 
and ABC [Fig. IHtL)] orientations. The different jumps 
and peaks in the DOS of Fig. HKc)-(d) are associated to 
the regions of the band dispersion marked by the hori- 
zontal red lines of Fig. [9l the energy of which depends 
on the values of the tight-binding parameters associated 
to inter-layer tunneling (71 and 73 in our case). In the 
two cases, we observe a splitting of the Van Hove peak, 
as seen in Fig. [8{e)-(f). Notice that when we have a high 
number of layers (e.g. above 10 layers), there is a weak 
effect on adding a new graphene sheet to the system, as 
it can be seen from the similar DOS between the 10- and 
the 50-layers cases of Fig. [8] 

In Fig. [10] we show the non-interacting (left pan- 
els) and the RPA (right panels) polarization function of 
MLG, for systems made of 3, 5 and 20 layers, and for 
ABA- and ABC-stacking. For the spectrum in the ab- 
sence of electron-electron interaction, as shown in Fig. 
[TOTa). (c) and (e), one does not observe any specific dif- 
ference in the two spectra apart from different intensities 
depending on the kind of stacking and on the number of 
layers considered for the calculation. On the other hand, 
the energy of the 7r-plasmon of ABC samples is rcdshiftcd 
with respect to the ABA-stacking. This can be see from 
the relative position of the peaks of — Im \(q,Lu) in Fig. 
fTOT b). (d) and (e). Also, notice that the separation be- 
tween the two peaks grows with the number of layers, and 
for a 20-layers system, the difference can be of the order 
of 1 eV, as it can be seen in Fig. ITOT f). In the follow- 
ing and unless we say the opposite, all the results will be 
calculated for the more commonly found ABA-stacking. 

For a more clear understanding about the evolution 
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FIG. 10. (Color online) Dynamical polarization and response 
function of ABA- and ABC-stacked multilayer graphene. The 
number of layers are 3 layers in (a,b), 5 layers in (c,d) and 
20 layers in (e,f). The size of each layer is 3600x3600 atoms 
for tri-layer (a,b); 3200x3200 atoms for 5-layer samples (c,d), 
and 1600x1600 atoms for 20-layer (e,f). 



of the particle-hole excitation spectrum with the number 
of layers, we plot in Fig. Illf a) the imaginary part of 
n(q, uj) for SLG and MLG of several number of layers, 
and compare the results to the polarization obtained us- 
ing the Dirac cone approximation. It is very important to 
notice that multi-layer graphene presents some spectral 
weight at low energies as compared to graphene, which 
can be seen from the finite contribution of ImII(q, uj) 
that appears to the left of the big peak of the graphene 
polarizability at uj — v-pq (~ leV for the used parame- 
ters), in terms of the Fermi velocity near the Dirac point, 
t>F = 3at/2. This is due to the low energy parabolic-like 
dispersion of bilayer and multilayer graphene, as com- 
pared to the linear dispersion of single layer graphene, 
and it can only be captured by considering the inter-layer 
hopping contribution to the kinetic Hamiltonian Eq. ([3]) . 
Furthermore, the spectrum presents a series of peaks for 
uj « v-pq, the number of which depends on the number 
of layers. This is due to the fact that as we increase the 
number of coupled graphene planes, the number of bands 
available for particle-hole excitations also grows leading 
to peaks at different energies for a given wave- vector 4^ 

The difference between SLG and MLG is also relevant 
in the low energy region of the dielectric function, as it 
can be seen in Fig. [LTT b). In fact, the ui — > limit 
of Re e(q, uj) calculated within the RPA grows with the 
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FIG. 11. (Color online) (a) ImII(q,a;) for SLG and MLG. The 
results for SLG obtained from the Kubo formula Eq. ^ are 
compared to those obtained from the Lindhard function Eq. 
(|lip . and to the Dirac cone approximation, (b) Re e(q, cj) for 
SLG and MLG, and comparison to the Dirac cone approxi- 
mation, for the same value of q as in (a). 



number of layers. Moreover, as we have discussed above, 
the zeroes of Re e(q, uj) signal the position of collective 
excitations in the system (plasmons). In Fig. ITTT b) 
we see that Re e(q, w), for the small wave- vector used, 
crosses for MLG, revealing the existence of a solution 
of Eq. dTBT) . but not so for SLG, as it was pointed out 
in Ref. |23| . However, we emphasize that the very exis- 
tence of solutions for the Re e(q, w) = equation for 
MLG docs not imply the existence of long-lived plas- 
mon modes. In fact, as we have already discussed in 
Sec. IIII1 these modes disperse within the continuum of 
particle-hole excitations [Im II(q, ui p i) ^ 0, where ui p i is 
the solution of Eq. (fT5]l]. so they will be Landau damped 
and will decay into electron-hole pairs with a damping 
given by Eq. (TTT)) . Furthermore, we remember that for a 
given wave-vector, the energy of the mode is controlled 
by the background dielectric constant k, as given by Eq. 
(f!?T]) . For the systems under consideration, n changes be- 
tween 1 (for SLG) and 2.4 (for graphite). The value of 
k, together with the form factor Eq. (TIT)]) that takes into 
account inter-layer Coulomb interaction, fix the position 
of the modes in each case. 

The effect of disorder in MLG is considered in Fig. 
IT2l where we show the polarization function of a 20-layer 
graphene system for different kinds of disorder. As in 
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FIG. 12. (Color online) Dynamical polarization and response 
function of ABA-stacked 20-layer graphene with disorders. 
The sample size of each layer is 1600 x 1600. 
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FIG. 13. (Color online) Loss function — Im l/e(q, u>) for SLG 
and MLG, which is proportional to the spectrum obtained by 
EEL experiments.— We have used the same q as in Fig. 1111 



the SLG, we find that disorder leads to a slight redshift 
of the peaks of the non-interacting spectrum [Fig. [P2Ta) 
and (b)], together with a smearing of the peaks at w ~ 
i>f9 and u ~ It. On the other hand, the interacting 
polarization function presents a reduction of the intensity 
of the plasmon peak due to disorder, as seen in Fig. UTTc) 
and (d), also in analogy with the SLG case. 



V. DISCUSSION AND COMPARISON TO 
EXPERIMENTAL RESULTS 

In this section we compute quantities which are di- 
rectly comparable to recent experimental results on SLG 
and MLG. We start by calculating the loss function 
— Im l/e(q, w), which is proportional to the spectrum 
measured by EELS. Our results, shown in Fig. 1131 are in 
good agreement with the experimental data of Ref. as 
in the experiments, we observe a redshift of the plasmon 



peak as one decrease the number of layers, as well as an 
increase of the intensity with the number of layers. No- 
tice that, due to finite size effects, there is an infra-red 
cutoff for the wave- vectors used in our calculations which 
prevents to reach the long wavelength limit. In Fig. Q2] 
we show the results for the smallest wave- vector avail- 
able, and we emphasize that the peaks will be further 
shifted to the left for smaller values of q. A further red- 
shift of the peaks would be obtained beyond RPA, as it 
has been reported for single- and bilayer graphene, where 
excitonic effects have been included^ 

We have also used our method to study the IXS ex- 
periments of Reed et alJz In Fig. [HJa)-(b) we plot the 
imaginary part of the non-interacting polarization func- 
tion for SLG and MLG, for two values of q similar to 
the ones used in Ref. As we have discussed in Sec. 
IIV1 inter-layer hopping leads to a finite contribution to 
the spectral weight in the low energy region of MLG as 
compared to the SLG spectrum. Notice that the number 
of peaks at this energy u « vpq scales with the number 
of accessible bands and therefore, with the number of 
layers. We emphasize that this effect is not included by 
the usually employed approximation of considering MLG 
as a series of single-layers of graphene, only coupled via 
direct Coulomb interaction^ Without the possibility of 
inter-layer hopping, the polarization function of graphene 
and graphite are, apart from some multiplicative factor, 
the same. As we have seen in Sec. IIV1 this simplification 
does not capture the low energy part of the spectrum, 
with some finite spectral weight due to low energy inter- 
band transitions between parabolic-like bands. 

At an energy of the order of u ~ 2t one observes the 
peak due to transitions between electrons from the Van 
Hove singularity of the occupied band to the singular- 
ity of the empty band. For SLG, the peak is split into 
two peaks if the wave-vector points in the T-K direction 
(as it is the case here), the separation of which increases 
with the modulus q = q x . However the amplitude of 
these peaks is highly suppressed from SLG to MLG. Fi- 
nally, one observes in Fig. mT b) that for higher values 
of q, as the one used here, there is a redshift of the peak 
of Im n(q,w) at the energy u ps VFq with respect to 
the Dirac cone approximation. This is due to trigonal 
warping effects, which are beyond the continuum approx- 
imation. Summarizing, we find two effects that lead to 
a global contribution to the polarization at low energies: 
one is the contribution to the spectral weight due to inter- 
layer hopping in MLG, and the other is the redshift of 
the peaks at w k v^q due to trigonal warping effects. 
Notice that, because we are studying the non-interacting 
polarization function, no excitonic effects are present in 
the results of Fig. H3Ta)-(b). 

Once the polarization function n(q, lo) is known, we 
compute the response function x(q, at the RPA level, 
as shown in Fig. [Htc)-(d). Again, we find a redshift 
of the position of the peaks as we decrease the number 
of layers. The different position of the peaks is due to 
the different contribution of inter-layer electron-electron 
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FIG. 14. (Color online) Non-interacting polarization function ImII(q, uj) [plots (a) and (b)] and RPA response function Imx(q, u>) 
[plots (c) and (d)] for SLG and MLG, for two different wave-vectors. The wave-vectors are chosen as in Ref. i. 
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FIG. 15. (Color online) Modulus of the screened fine structure 
constant |q*|, calculated from Eq. (|23|l . for SLG and MLG 
of a different number of layers. The inset is a zoom for the 
more experimentally relevant uj — > region of the spectrum 
(see text). \a*\ ~ 0.6 for SLG, whereas this value is highly 
reduced for MLG: |a*| ~ 0.3 for a 20-layers sample in our 
numerical calculation, which has a behavior very similar to 
graphite. 



interaction for each well as to the different value 

of k as given by Eq. (|2"Tj) . Our results agree reasonably 
well with those of Ref. [5J 

Finally, we calculate the renormalization of the fine 
structure constant a = e 2 /v-p due to dynamic screening 



associated to the inter-band transitions from the valence 
band. For this, and in analogy with Ref. we define 



a*(q,w) 



e(q,w) 



(23) 



The results for the modulus |a*| for SLG and MLG are 
shown in Fig. [TS] In this plot we have used the value of 
the Fermi velocity valid near the Dirac point, i.e., vf = 
3at/2. Therefore, we emphasize that these results should 
be reliable only at low energies. At u> — > and for the 
smallest wave-vector we can access (q = 0.2a _1 ), RPA 
predicts |a*| rj 0.6 for SLG, which is considerably higher 
than the value estimated in Ref. [f| |a*| sa 0.15. However, 
the results that we obtain for MLG are much closer to this 
value: we find that |a*| ~ 0.3 for graphite, only slightly 
higher (a factor of 2) than the experimental results of 
Ref. i, which are actually obtained from graphite. 



VI. CONCLUSIONS 

In conclusion, we have studied the excitation spectrum 
of single- and multi-layer graphene using a full 7r-band 
tight-binding model in the random phase approximation. 
We have found that, for MLG, the consideration of inter- 
layer hopping is very important to properly capture the 
low energy region {uj ~ v-pq) of the spectrum. This, to- 
gether with trigonal warping effects, lead to a finite con- 
tribution to the spectral weight at low energies as well 
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as a redshift of the peaks with respect to the Dirac cone 
approximation. We have also studied the high energy 
plasmons which are present in the spectrum of SLG and 
MLG at an energy of the order of cj ~ 2t ps 6eV and 
which are associated to the enhanced DOS at the Van 
Hove singularities of the 7r-bands. The energy of the 
7r-plasmon depend also on the orientation between adja- 
cent layers, and we find that, for a given wave- vector, the 
energy of the mode for ABC-stacked MLG is redshiftcd 
with respect to the corresponding energy of ABA order- 
ing. This difference is higher as we increase the number 
of graphene layers of the system. 

The effect of disorder has been considered by the in- 
clusion of a random on-site potential and by a renormal- 
ization of the nearest neighbor hopping. Both kinds of 
disorder lead to a redshift of the ui rts vpq and u) « 2t 
peaks of the non-interacting excitation spectrum and to 
a smearing of the Van Hove singularities. The position 
of the 7r-plasmons is unaffected by disorder, although the 
height of the absorption peaks is reduced as compared to 
the clean limit. 

Finally, we have compared our results to some re- 
cent experiments. Our calculations for the loss function 



Im l/e(q, to) show a redshift of the SLG mode with re- 
spect to graphite, and compare reasonably well with ex- 
perimental EELS dataj&^ Furthermore, we also obtain 
good agreement with the IXS results for the response 
function obtained in Ref. [^. We obtain a static dielectric 
function which grows with the number of layers of the 
system. In the long wavelength and lu — > limit, the 
dynamically screened fine structure constant is found to 
be highly reduced from graphene to graphite. The value 
that we find for a MLG in the RPA, without considering 
any excitonic effects, is about two times larger than the 
one estimated in Ref. [H for graphene. More accurate re- 
sults could be obtained going beyond single-band RPA^ 
which is beyond the scope of this work. 
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